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FOURIER  TRANSFORM  RECONSTRUCTION  FROM  INEXACT  DATA 


INTRODUCTION 

The  central  problem  in  spectral  estimation,  the  reconstruction  of  a  Fourier  transform  from 
sampled  data,  is,  to  paraphrase  Parzen  [1] ,  essentially  the  problem  of  how  best  to  approximate  a 
function  from  N  of  its  Fourier  coefficients.  Emphasis  on  the  approximation-theoretic  aspect  of  the 
problem  focuses  attention  on  the  algebraic  form  of  the  estimating  functions  and  on  the  error  crite¬ 
rion  by  which  the  approximation  is  judged.  In  [2]  Fourier  transform  reconstructions  were  ob¬ 
tained  from  noise-free  data  samples  by  using  best  linear  approximations  in  suitably  chosen  Hilbert 
spaces.  These  reconstructions  included,  as  special  cases,  the  estimates  obtained  by  the  bandlimited 
extrapolation  procedures  of  Cadzow  [3] ,  Papoulis  [4] ,  Kolba  and  Parks  [5] ,  and  others.  The  same 
approach  was  used  in  [6]  to  derive  algorithms  for  tomographic  reconstruction.  In  this  report  we 
consider  the  problem  of  Fourier  transform  estimation  from  data  corrupted  by  additive  noise  and 
we  analyze  the  reconstructions  obtained  in  [2]  as  statistical  estimators. 

We  begin  with  a  summary  of  the  techniques  developed  in  [2]  for  the  noise-free  case.  These 
methods  are  then  extended  to  the  case  of  noisy  data,  and  bias  and  variance  formulae  are  obtained. 
The  special  case  of  sinusoids  in  noise  is  considered,  and  the  amplitude  estimates  obtained  for  this 
case  are  shown  to  be  superior  to  the  best  linear  unbiased  estimates. 

The  general  problem  to  be  considered  here  is  the  reconstruction  of 

F(w)=  f  f(t)  eiwt  dt ,  (1) 


given  the  values  at  t  =  , ...,  tN  of 


*(*)  =  W)  +  n(t).  (2) 

Specific  assumptions  about  n(f)»  the  unwanted  or  noise  component,  will  be  made  as  needed.  By  the 
noise-free  case  we  will  mean  x(t)  =  f(t). 


APPROXIMATING  A  FOURIER  TRANSFORM  FROM  NOISELESS  DATA 

The  problem  of  estimating  F(w),  given  f{t± ), ...,  f(tN),  is  often  referred  to  as  the  quadrature 
problem;  one  must,  in  effect,  perform  the  integration  in  Eq.  (1)  in  some  sense.  Unless  more  than 
just  the  numerical  data  are  known  about  F(w)  (or  equivalently,  about  f(t)),  any  function  F(w)  such 
that 


A*.) 


-/ 


F(w)e~,wtn  dw/2n  ,  n  =  l,...N, 


(3) 
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is  a  possible  solution.  The  main  idea  of  [2]  is  that  one  must  use  whatever  prior  knowledge  one  has 
to  eliminate  solutions  as  unsuitable  and  to  select  a  single  F(w )  that  stands  the  best  chance,  as  far 
as  one’s  prior  information  is  concerned,  of  being  the  true  F(w).  An  effective  way  of  doing  this  is 
by  the  use  of  weighted  L2 -spaces  and  best  linear  approximation. 

Let  F  (w)  be  a  prior  estimate  of  F(w)  and  let  P(w)  >  0  be  a  weighting  function.  We  take  as 
the  estimate  of  F(w)  that  function  F(w)  that  minimizes  the  weighted  error 

f  [)F0(u>)-  F(w)\2IP(w)]dw  ,  (4) 

subject  to  the  data  constraints  of  Eq.  (3).  The  solution  is  easily  shown  to  be 


N 

F(w)  =  FQ(w)  +  P(w)  £  *meiwtm  .  (5) 

m“  1 

where  the  coefficients  zlt zN^ are  determined  by  Eq.  (3).  In  addition  to  being  closest  to  Fa  in  the 
sense  of  Eq.  (4),  this  estimate,  F(w),  is  the  best  estimate  of  F(w)  having  the  algebraic  form  of 
Eq.  (5)  in  the  sense  that  the  error 

oo  JV 

f  [>(«;)-  jF0(u>)  +  P(u>)  f^ameiwtm\\2/P(w)]dw,  (6) 

J-aa  m  =  l 

is  a  minimum  when  am  =  zm ,  m  =  1, ....  N.  Although  the  use  of  a  nonzero  Fa(w)  provides  added 
flexibility  when  such  a  prior  estimate  is  available,  we  will  take  F0(w)  -  0  in  what  follows.  The  esti¬ 
mator,  Eq.  (5)  with  Fa(w)  =  0,  was  considered  in  [2] ,  where  it  was  referred  to  as  the  PDFT 
estimator. 

If  the  data  are  evenly  spaced  d  units  apart  and P(w)  =  1  for  lu>l<7r/d,P(u>)  =  0  otherwise,  then 
Eq.  (5)  reduces  to  the  DFT  of  the  data.  If  we  should  happen  to  know  that  F(w)  =  0  for  \w  |  >  o, 
where  0  <  a  <  ir/d,  and  we  take  P(w)  =  1,  \w  \  <  a,  P(w)  =  0  otherwise,  we  obtain  the  estimate 


Fa(«,)  = 


N 

E  zme 

m  =1 
0 


iwmd 


Iwl  <  a 
|u;|  >  a 


(7) 


where  zl ,  z2, ...,  zN  are  determined  by  Eq.  (3).  As  was  shown  in  [7]  this  estimate  is  the  closed-form 
equivalent  of  those  obtained  through  extrapolation  algorithms  in  [3] ,  [4] ,  [5]  and  elsewhere.  The 
estimator  in  Eq.  (7)  is  quite  unstable,  based  as  it  is  on  the  assumption  that  the  data  correspond  to  a 
signal  that  is  precisely  o-band-limited.  It  was  in  an  attempt  to  reduce  the  sensitivity  exhibited  by 
Fa(«>)  to  out-of-band  components  that  we  were  led  to  consider  the  estimators  in  Eq.  (5).  We  found 
that  by  taking  P(w)  =  1  for  |  w  |  <  a,  and  P(w)  =  e  >  0  for  a  <  |  w  i  <  ir/d,  where  e  is  a  small  positive 
number,  the  resulting  PDFT  estimator  provided  a  marked  improvement  over  F0(w).  In  effect,  we 
make  P(w)  share  with  the  true  F(w)  broad  features,  such  as  we  know  them  a  priori.  Relative  energy 
concentrations  and  known  component  distributions  such  as  radar  clutter,  noise,  and  even  delta  func¬ 
tions,  can  be  incorporated  in  P(w).  As  we  begin  to  distinguish  components  from  one  another  and  to 
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take  them  into  account  in  P(w),  we  pass  into  what  may  properly  be  called  the  noisy  case  and  this 
is  the  subject  of  the  next  section. 


STATISTICAL  ESTIMATION  OF  FOURIER  TRANSFORMS 

Suppose  now  that  x(t)  =  f(t)  +  n(t)  and  that  ), ....  x(tN)  are  given.  We  shall  obtain  a  PDFT 
estimate  of  X(u>),  the  Fourier  transform  of  x(t),  and  then  use  it  to  derive  an  estimate  of  F{w).  We 
select  weighting  functions,  PAw)  >  0  and  Pn  (w)  >  0  (which  correspond  to  f(t)  and  n(t)  respec¬ 
tively),  take  P(w)  =  Pf(w)  +  Pn(w),  and  obtain 


N 

X(w)  =  P(w)  £  2meiwtm  ,  (8) 

m  =  l 

as  the  estimate  for  X(u;).  As  a  reasonable  choice  for  the  estimate  of  F(w)  we  choose 

F(w)  =  [pf(w)/P(w)]x{w) , 


F(w)  =  Pf(w)  (9) 

m  =  1 

where  the  zx , ...,  zN  are  chosen  to  make  Eq.  (8)  consistent  with  the  data.  In  the  special  case  in 
which  Pf(w)Pn(w)  -  0  (disjoint  spectra),  it  is  clear  that  X(u>)  estimates  both  F(w)  and  N(w),  the 
Fourier  transform  of  n(t).  An  example  of  this  case  was  considered  in  [2]  where  Doppler  radar 
samples  were  analyzed  to  detect  targets  known  to  be  outside  the  clutter  spectrum. 

To  compute  the  bias  and  variance  of  F(w)  treated  as  a  statistical  estimator  of  F(w),  it  is  nec¬ 
essary  to  make  further  assumptions  about  n(t).  Before  doing  that,  however,  it  is  worth  noting  that 
up  to  now  n(f)  has  been  any  additive  component  that  we  desire  to  remove  by  filtering.  If  we  take 


OO 

t)=  f  F(w)e~  iwtdw/2ir , 

J—  OO 


N 


=  E  V/<f-  fm)’ 

m  *  1 


(10) 


where  Pf(t)  is  the  inverse  Fourier  transform  of  Pf{w),  then  we  have  filtered  x(t).  Indeed,  by  setting 


N 


ft*,,)-  E  zmPfVn'  n  =  N, 

m  •  1 


(ID 


and  then  usingAthese  values  as  noiseless  data  in  Eq.  (5)  (with  F0(w)  =  0,P(u;)  =  Pf(w)),  we  reproduce 
the  estimate,  F(w),  of  Eq.  (9).  The  PDFT  can  in  this  way  be  used  to  introduce  prior  information 
into  linear  filtering. 

From  now  on  it  will  be  assumed  that  n(t)  is  one  realization  of  a  mean-zero  weakly  stationary 
random  process,  with  power  spectral  density  function,  Rn(w),  and  autocorrelation  function,  r„(f)- 
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With  P(u>)  =  Pf(w)  +  Pn(w)  and  p(t)  the  inverse  Fourier  transform  of  P(w),  let  G  be  the 
N  XN  matrix 


G  =  [p(tn-  fm)]  , 

and  let  C  be  its  inverse 

C-lCn<m)=G-K 


Then  we  can  write 


F(w) 


N 


m  =  1 


(12) 


(13) 


(14) 


with 


M 

Hm(w)-Pf(w)YJCnmeiu3ttn,  m  =  1 . N.  (15) 

n»l 

It  will  be  convenient  to  define,  for  any  vector  v  =  (i^ , ....  vN), 

M 

T(v,w)=  Y,  Hm(w)vm  .  (16) 

m  "  1 

The  expected  value  of  F(w)  is  then 

E{F(w))  =  T(f,w),  (17) 

for 


so  that  the  bias  is 


F{w)  -  E(F{w))  =  F(w)  -  T(f,w) .  (18) 

Note  that  T(f,w)  is  the  estimate  of  F(w)  that  would  be  obtained  if  the  data  vector  were  the  noise¬ 
less  vector  f.  One  source  of  bias  is  obviously  traceable  to  deficiencies  in  the  prior  information  even 
in  the  noiseless  case  (the  quadrature  problem).  Another  source  of  bias  stems  from  the  decision  to 
design  the  estimator  F{w)  to  operate  in  a  noisy  context  and  to  use  a  prior  estimate  of  the  signal- 
to-noise  ratio  of  roughly  Pf(w)/Pn(w).  The  bias  in  Eq.  (18)  depends  solely  on  what  would  happen 
if  the  noise  mysteriously  vanished.  We  purposely  introduce  bias  so  that  the  variance  can  be  reduced. 

The  variance  of  F(w)  is  given  by 


N  N  _ 

var[F(u>)]=£  £  HmW  Hn^P^m  ~  ’ 

m  - 1  n  *  1 


(19) 
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or 

vai[F(w)}=  f  \T(e(d),w)\2Rn(d)d6/2ir ,  (20) 

oo 

where  e(d)  =  (exp  (-it18), ...,  exp (-itN0)).  Note  that  T{e(6),w)  is  the  value  of  the  estimator  F(w) 
based  on  a  data  vector  corresponding  to  samples  of  a  pure  sinusoid  at  frequency  6 .  To  the  extent 
that  Pf(w)  effectively  describes  F(w),  it  is  to  be  expected  that  most  of  the  e(6)  data  will  be  attrib¬ 
uted  to  noise,  especially  if  Rn(d)  is  large  and  this  is  reflected  in  Pn(6).  The  spillover  of  noise  at  6 
into  signal  estimated  energy  at  w  should  be  small. 

Next  we  consider  the  special  case  of  sinusoids  in  noise. 


SINUSOIDS  IN  NOISE 

Let  the  information-bearing  component  of  the  received  signal  consist  of  K  sinusoids 

K 

m-  E  (21> 

fc= l 

and,  as  before,  let  the  given  data  be  x(fn ),  n  -  1, ....  N,  where  x(f)  =  f(t)  +  n(t).  The  problem  posed 
is  to  estimate  the  number  of  sinusoids  present,  K,  their  frequencies,  w1 , ...,  wK,  and  their  complex 
amplitudes,  bj , ...,  bK. 

No  single  weighting  function  is  suitable  for  determining  both  the  locations,  wk,  and  the  ampli¬ 
tudes,  bk,  of  the  sinusoids.  Before  the  number  and  location  of  the  sinusoids  are  known,  Pf(w) 
should  be  chosen  to  be  flat,  corresponding  to  our  knowledge  that  almost  all  of  F(w)  is  zero.  The 
estimate,  F(w),  thus  obtained  then  provides  information  about  where  spectral  peaks  are  not  to  be 
found.  Where  peaks  are  located  we  can  expect,  using  this  flat  Pf(w),  to  see  only  moderate  indica¬ 
tions  in  F(w),  because  of  the  algebraic  limitations  of  the  estimating  function  (9).  After  a  good 
estimate  of  the  peak  locations  has  been  obtained,  that  information  can  be  incorporated  in  the 
estimator  by  choosing  a  new  weight,  Pf(w),  which  is  larger  in  the  neighborhood  of  likely  sinusoidal 
peaks.  In  the  limiting  case,  when  the  values  of  wk  are  known  exactly,  a  weighting  function  contain¬ 
ing  delta  functions  can  be  used  to  estimate  the  complex  amplitudes,  bk . 

To  consider  this  limiting  case  in  more  detail,  suppose  the  sinusoid  frequencies,  w1 , ....  wk, 
have  been  determined  and  Pf(w)  has  been  set  equal  to 

K 

Pf(w)  =  2tt  23  8(w-wk) ,  (22) 

fc-l 

so  that 

K 

•  (23) 

ft  =  i 
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The  noise  component  is  associated  with  a  density  Pn(w),  whose  inverse  transform,  p, 


where 


Pn(0  =  “*„(f), 

a=Pn(0)>0  . 


The  matrix  G  in  Eq.  (12)  becomes 

=  Ip  fit j  ~  tm)  ]  +  «[*n(fy-  tm)  ]  . 

Letting  J be  the  N  X  K  matrix  J  =  [exp  (~iwktm )]  we  have,  for  Gn  =  [gn(t:  -  fm)], 

G  =  J/"  +  oc  Gn  . 

Using  Eq.  (22)  in  Eq.  (9)  we  can  write  F(w)  as 

K  N 

£(£ 
fc  =  l  m  =  1 

and,  thus,  from  Eq.  (21)  the  estimate  for  bk  is  given  by 

N 


F(w)  =  2ir  £  zme,Wf*tm)  &(,w~  wk) , 


bk  "  £  zme'w^tm  , 


m  -1 


where 

and  where 


Letting  6T  *  (6lf ...,  bK),  we  obtain 


and 


zr  =  (z1,...,*JV)  =  xrG  1  , 
X  —  (x(tj  ),  ...,  x(fjy))  . 

b  =  J*z, 


b  =  j*G~1x 
-  +  OC  Gn)~1x 

=j*(g;1jj*  +  oc/)-ig-ix, 

where  /  is  the  identity  matrix.  A  simple  calculation  shows  that 

+  <x/)-l  =  («/  +  J*G“1«f)~1J*  , 


(t),  we  write  as 


(24) 

(25) 

(26) 

(27) 

(28) 

(29) 

(30) 
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and  so  we  obtain, 


b=(<xi  +  j*G~1J)'1J*G~n1x.  (31) 

It  follows  that  as  a  converges  to  zero,  our  estimating  vector  b  converges  to 

bo  =  (J*G-nlJ)-'j*G-n1x,  (32) 

which  is  easily  recognized  as  the  best  (minimum  variance)  linear  unbiased  estimate  (BLUE)  of  the 
coefficients  , ...,  bK .  By  taking  a  to  zero,  we  gradually  eliminate  the  noise  component  of  the 
prior  P(w).  As  we  saw  earlier,  this  noise  component  is  the  source  of  bias,  and  so,  in  the  limit  when 
a  goes  to  zero  we  obtain  an  unbiased  estimate.  The  estimate  in  Eq.  (28)  is  biased  to  account  for  the 
assumed  presence  of  a  noise  component  at  a  level  of  relative  to  the  signal  energy  level  described 
by  Pf(w).  If  the  prior  assessment  of  SNR  (signal-to-noise  ratio)  is  accurate,  it  is  to  be  expected  that 
Eq.  (28)  will  provide  a  better  estimate  than  the  BLUE. 

The  inversion  of  the  matrix  G  in  Eq.  (24)  is  simplified  by  observing  that  JJ*  is  a  sum  of  outer 
products  (column  matrix  times  row  matrix).  An  efficient  scheme  based  on  the  identity 


( A  +  xxT)~ 1 


..i  (A-1x)(*tA"1) 

A - 

1  +  xtAx 


for  column  vector  x,  has  been  used  successfully  in  simulations  we  have  run. 


SUMMARY 

We  have  described  methods  for  the  linear  reconstruction  of  Fourier  transforms  from  noise¬ 
less  data  and  have  extended  them  to  the  case  of  noisy  samples.  Bias  and  variance  of  these  estima¬ 
tors  in  the  presence  of  weakly  stationary  random  noise  were  calculated  and  the  special  case  of 
sinusoids  in  noise  was  considered.  We  found  it  useful  to  view  this  latter  problem  as  one  of  fre¬ 
quency  estimation  followed  by  amplitude  estimation.  The  resulting  biased  amplitude  estimates  were 
compared  to  their  limiting  unbiased  values,  which  were  shown  to  correspond  to  the  BLUE. 
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